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The finite size corrections in exact diagonalization studies 
of the one- and two dimensional Hubbard model can be re- 
duced systematically by a grand-canonical integration over 
boundary conditions. We present results for ground-state 
properties of the 2D Hubbard model and an evaluation of 
the specific heat for the ID and 2D Hubbard model. We find 
the reduction of the finite size corrections to be substantial, 
especially in 2D. 
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INTRODUCTION 

The exact diagonalization of the Hubbard model (or 
the t — J model) on finite clusters has become one of the 
promineiit techniques in the study of correlated electron 
systemsEl. This technique has the advantage of treating 
the correlations on finite clusters without any approxima- 
tion. Since the Hilbert space grows exponentially with 
cluster size one is limited, for the Hubbard model, to 
clusters with up TVs — 10 sites. By finite size corrections 
one means the difference between the results obtained for 
the finite cluster and its (unkown) value in the thermo- 
dynamic limit, Ng = oo. 

The control of the finite size corrections poses a great 
challenge to the exact diagonalization technique, espe- 
cially in the case of two dimensions (2D). In ID, it is 
a standard procedure to plot the results systematically 
as a function of chain length, Ng = 6, 8, 10, 12, . . . sites. 
In 2D, on the other side, where most calculations have 
been done up to now with either periodic or antiperi- 
odic boundary conditions, it is difficult to estimate the 
finite size corrections by direct comparison of clusters of 
different sizes, for two reasons: 

1. The finite size corrections are often non-monotonic 
as a function of Ng , due to strongly varying cluster 
geometries. 

2. The nominal density of particles, n — Nf,/Ns does 
normally not coincide for clusters of different sizes 
as the allowed number of particles Ng = 1,2,3,... 
is an integer. 

Both of these two difficulties can be circumvented by a 
method introduced recentlya, the integration over bound- 
ary conditions (IBC). The IBC circumvents the first of 
above difficulties by performing a grand-canonical inte- 
gration over boundary conditions. It has been shown. 



that this procedure removes all those finite size correc- 
tions which are causedp.by the special geometry of the 
Fermi-sea of the clusteiu. The grand- canonical approach 
then leads to the possibility of directly comparing the 
results for different Ng with the same density n, solving 
also the second of above difficulties. 

In a previous publicationH the IBC had been intro- 
duced and some results for the Hubbard model on very 
small clusters had been presented. Here we want to study 
in detail how the finite size corrections can be estimated 
systematically within the IBC both at zero and at finite 
temperatures. For this purpose we will concentrate on 
the Hubbard model, the same study could be performed, 
in principle, for the t ^ J model or any other cluster 
Hamiltonian. 



INTEGRATION OVER BOUNDARY 
CONDITIONS 

We consider the Hubbard Hamiltonian on a cluster 
with sites, 



<i,j>,(y 



(1) 



where the symbol < i,j > denotes pairs of n.n. sites 
and where the c\ ^ and the c^ ^ are the electron cre- 
ation/destruction operators of spin a =1,1 respectively. 
The hopping amplitudes. 
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depend, in general, on phases onj, related to the bound- 
ary condition. In the following all energies will be mea- 
sured in unities of t. We will now give an account of the 
IBC for ID, the generalization to 2D is then straightfor- 
ward. 

In ID we can choose the a^j — a/Ng, where a G [0, 27r[ 
is the boundary condition. Periodic and antiperiodic 
boundary conditions correspond to a = and a = tt 
respectively. The IBC technique needs the exact diago- 
nalization of the cluster for any particle number Ng = 
0,1,2,..., 2Ns and any a S [0, 27r[. For any given a one 
then calculates the free energy 



F,(/i) = -TkB ln{ 
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where /3 = l/^ksT) is the inverse temperature, /i the 
chemical potential, Nh ~ 1 the dimension of the Hilbert 
space and the Ek{a, Ne) the eigenenergies. {Nh depends 
on both on Ng and Ne). The total free energy of the 
cluster and the particle density n S [0, 2] are then given 
by the average 

Fif^)- (4) 

Jo 

and by n = The free energy Eq. (jj) ob- 

tained with the IBC technique is constructed in such a 
way that the firHte size corrections vanish identically in 
the limit U — 00. One can understand this result most 
easily in momentum space, noting that the Fermi -sea, as 
obtained by the IBC, has no finite size corrections for 
U = 0. This property of the IBC holds in any dimension 
and is especially important in 2D. The particle density 
n can furthermore be tuned to any value, by choosing 
an appropiate fi, allowing to directly compare data from 
clusters with different Ns. These properties distingish the 
IBC technique from other exact diagonaliaation studies 
employing a range of boundary conditionsa. 

The integration over boundary conditions, appearing 
in Eq. (^, is replaced in actual calculations by a fi- 
nite sum over Na equal-distant boundary conditions. In- 
deed we will find further below that only a few Na are 
needed in order to already obtain most of the reduction 
in the finite size corrections. Ground-state properties 
can be calculated, as usual, by the Lanczos technique. 
With Eo{a, N^) being the ground-state energy for a given 
a and TVg one finds the estimate eo(/i) for the grand- 
canonical ground-state energy per site to be 

s a ^ 

A Legendre transformation then yields the canonical 
ground-state energy of the cluster, eo{n,Ns), as a func- 
tion of density. In Fig. |l| we have plotted results obtained 
for the 2D Hubbard model with U = At. The results for 
Na = 1, corresponding to periodic boundary conditions 
only, differ typically by 10% — 20% in between clusters 
with Ns ^ 8 and N^ = 10 sites. -Also shown in Fig. |l| 
are the data obtained for Na — 250. Here the data differ 
only by 1% — 2% in between Ng ~ 8 and Ng — 10, an im- 
provement by one order of magnitude. For comparision 
we have plotted also data obtained by projection Monte 
Carlo for a 10 x 10 system with periodic boundary con- 
dition by Furukawa and ImadaB. 

SPECIFIC HEAT 

Eq. (||) for the free energy needs the knowledge of all 
eigenstates and can therefore be used only for very small 
sytems, typically up to A^s = 6 for the Hubbard model. 
Several approximative methods have been developed for 



the numerical evaluation of theriaadynamic properties of 
sytems with large Hilbert spacesErO. Here we use the ker- 
nel polynomial M,pproximation (KPA) developed by Sil- 
ver and RoederQ, since it has a very good error control. 
One starts by scaling the Hamiltonian so that the mag- 
nitude of all eigenenergies is less than unity. The grand- 
canonical partition function, Za{^i) = exp[— /3_F'c((^)], is 
then expressed as an integral over the density of states, 
Daiuo.Ne) = llNHY.k=o^ 5{u ~ Ek{a,Ne)), 

2Ns 1 

Zait-i) = ^^e^^ / dujDa{iu,Ne)e-'''^. (6) 

The density of states can be expanded in a set of or- 
thogonal polyomials, here we take Legendre polyomials, 

Da{uj,Ne)^ 2 Mc^,Ne)Pi{uj), (7) 

with Njnom — *■ oo. Using the orthogonality relations 
for Legendre polynomials one can express the coefficients 
ai{a,Nf.) in term of the eigenstates of the Hamiltonian, 

H\Ek{a,N,)) ^ Ek{a,Ne)\Eu{a,Ne)), via 

az(a,Ae) = — — {Ek{a,Ne)\Pi{H)\Ek{a,Ne)). 

^ k=Q 

(8) 

In practice one uses a finite Nmom < oo. A finite values 
of Nrnom Icads to the well known Gibbs oscillations in 
Da{uj, Nf.). The Gibbs oscillations can be smoothed out 
by the replacement a; — > a; gi{Nmom) in (^. The optimal 
functional form of the smoothing functions gi{Nmom), 
with go{Nmom) = 1 and 5w„o„+i(A^mom|)| = 0, has 
been studied intensively in the literaturelD. We use 

giiNmom) = {sin{z)/zY, with Z = lTr/{Nmom + !)• 

The number Nmom of polynomials necessary for an ac- 
curate representation of the density of states in (|^) in- 
creases with decreasing temperature, T. We have evalu- 
ated the specific heat cy — (3'^ < {H— < H >)^ > for 
a 6-site Hubbard chain with periodic boundary condi- 
tions exactly, via Eq. (^), and via the kernel polynomial 
approximation with various Nmom. We have plotted in 
Fig. I the results for (a) t/ = 2 and (b) U = 16. F 
The KPA becomes asymptotically exact for large tem- 
peratures and any values of Nmom,. A numerical accu- 
rate approximation to the specific heat for temperatures 
down to r ^ 0.25t may be obtained with Nmom ^ 10'^ 
and Nmom ~ lO'*, for [/ = 2 and U — \& respectively. 
The large value, Nmom ^ iC, needed for big C/'s is the 
reason that all data presented further below will be for 
U = 2. 

Formula (||) is clearly not useful for larger clusters. 
It has been observedQ that one can systematically ap- 
proximate the trace occuring in the RHS side of (||) by 
sampling over Nr random states |r). 
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m{a,N,)^^Y.^r\Pi{H)\r)- (9) 

^ r— 1 

It is possible to evaluate recursiV|ely, by a procedure 
very similar to the Lanczos techniqueQ. Theeprors intro- 
duced by random averaging vanish like 1 / v^Ayj. We have 
performed extensive tests of the dependence of the data 
on Nj.] we will present further below only data which are 
unaffected by finite- -/V^ effects. 

RESULTS 

In Fig. H we have plotted the results for the spe- 
cific heat per site for various ID chains of length Ng = 
4, 6, 8, 10, both for the case of periodic boundary condi- 
tions only, Na = 1, and for Na = 10. Let us first note 
that the specific-heat curves for clusters of different sizes 
have eventually all to coincide at large enough tempera- 
tures, due to the grand-canonical formulation. At large 
enough temperatures the leading terms contributing to 
the specific heat of a cluster of a given size are identical 
to the leading terms of a high-temperature expansions. 
One can consequently define a T*{Ns, e) as the tempera- 
ture above which the finite size corrections to the specific 
heat are smaller than a given accuracy e, say e 1%. As 
the specific heat for iVg = oo is generally not known we 
took as a practical estimate for T* {Ns , e) the criterium 
that \cy (T, N,)-Cy {T, iV, - 2) | < e for aU T > T* {Ns , e) . 
We found that T* {Ng, e) is roughly inversely proportional 
to Ns, for the ID data presented in Fig. ^ and that the 
integration over boundary conditions reduces T* (Ng , e) 
by factors of 2 — 3, as it is evident from a comparison of 
Fig. |(a) and Fig. |(b). 

The IBC improves the cluster estimates of the specific 
heat in a second way, besides the reduction of T*{Ns, e). 
One finds empirically that the specific heat becomes lin- 
ear at low temperatures in the limit Na — *■ oo. This prop- 
erty of the low-T specific heat is present for all Ns and is a 
consequence of the absence of those finite size corrections 
which are caused by the geoiaetry of the Fermi-surface 
of the cluster, within the IBCcl. For finite Na there will 
be in general an intermediate temperature range where 
the specific heat is roughly linear, like it is the case for 
A^'^ = 10 in Fig. |(b) for temperatures O.lt - T - OM. 
This feature of the specific heat data obtained with the 
IBC may be used, e.g., to estimate the inverse effective 
mass. 

In Fig. ^ we present the specific heat per site, as a 
function of temperature, for 2D clusters with Ns = 8, 10, 
U = 2t,n = 0.8 and Na = l and Na = 16. The improve- 
ment obtained with the IBC is even more pronounced 
than in ID. The data for periodic boundary conditions is 
affected so much by finite size errors that is does not allow 
for a reliable estimate of the specific heat for T < 2t. For 
Na = 16 the finite size corrections are, on the other side, 
practically absent for T > T*{Ns = 10, e - 1%) ~ 0.72i. 



In conclusion we have shown that the integration over 
boundary conditions (IBC) can be used to reduce sub- 
stantially the finite size corrections occuring in exact di- 
agonalization studies. The data obtained by the IBC 
show typically a smooth behaviour as a function of clus- 
ter size, even in 2D where the cluster geometry may vary 
widely from cluster to cluster. This smooth behaviour 
often allows for quantitative estimates of the finite size 
corrections of the data obtained by the IBC, even in 
2D. One obtains furthermore certain qualitative improve- 
ments with the IBC, like a linear specific heat for the 
Hubbard model on one- and two dimensional clusters. 
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FIG. 1. The ground-state energy, eo, per site as a func- 
tion of particle densitiy n. Plotted are the results for the 2D 
Hubbard model and U — 4t. The filled circles and squares 
denote the results for clusters with periodic boundary condi- 
tions only, Na = 1, for clusters with Ns = 8 and Ns — 10 
sites respectively. The long-dashed and the continuous line 
are the respective results for Na = 25. Also shown (stars) is 
data obtained by projection Monte Carlo for a 10 x 10 system 
by Furukawa and ImadaEI. 

FIG. 2. The specific heat, cy, per site, as a function of 
temperature, of a 6-site chain with periodic boundary con- 
dition, n = 0.8 and (a) U ^ 2t and (b) U = 16t. Plotted 
are the exact results (solid line) and the results of the kernel 
polyomial approximaiton with (a) Nmom ~ 10, 100, 1000 and 
with (b) Nmom = 100, 1000, 10000. 
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FIG. 3. The specific heat, cv, per site as a function of 
temperature, for ID chains with Ns = 4, 6, 8, 10, U = 2t, 
n = 0.8 and (a) Na = 1 and (b) Nc = 10. For iV^ = 4, 6 the 
data are exact, for = 8, 10 the KPA has been used with 
Nmom = 1000 and Nr = 100, 10 for = 8, 10 respectively. 



FIG. 4. The specific heat, cv, per site as a function of 
temperature, for 2D clusters with — 8, 10, U = 2t, n = 0.8. 
Plotted are the data for Na = 1 (dotted/dashed line) and 
Net = 16 (dashed-dotted/solid line). The KPA has been used 
with Nmom = 1000 and Nr = 10, 8 for A^q, = 1, 16 respectively. 
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